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ABSTRACT 

We study the spatial distribution of Galactic metal-free stars by combining an extremely high- 
resolution (7.8 x 10 5 Mq per particle) Cold Dark Matter N-body simulation of the Milky- Way with 
a semi-analytic model of metal enrichment. This approach allows us to resolve halos with virial 
temperatures down to the 10 4 K atomic cooling limit, and it is sufficiently flexible to make a number 
of robust conclusions, despite the extremely uncertain properties of the first stars. Galactic metal-free 
stars are formed over a large redshift range, which peaks at z ~ 10, but continues down to z ~ 5, 
contributing stars at wide range of Galactocentric radii. Stars containing only metals generated by 
primordial stars are similarly widespread. Neither changing the efficiency of metal dispersal by two 
orders of magnitude, nor drastically changing the approximations in our semi-analytical model can 
affect these result. Thus, if they have sufficiently long lifetimes, a significant number of stars formed in 
initially primordial star clusters should be found in the nearby Galactic halo regardless of the specifics 
of metal-free star formation. Observations of metal abundances in Galactic halo stars should be taken 
as directly constraining the properties of primordial stars, and the lack of metal-free halo stars today 
should be taken as strongly suggesting a 0.8Mq lower limit on the primordial initial mass function. 
Subject headings: Galaxy:formation - Galaxy:evolution - stars:abundances - cosmology: theory 



1. INTRODUCTION 

Metals are made in stars, but stars need not to have 
metals. While the early universe was an efficient source 
of helium, primordial nucleosynthesis was halted before it 
was able to produce elements heavier than lithium. Thus 
the first stellar generation of Population III (PopIII) 
stars condensed from a gas that was free of carbon, sil- 
icon, iron, or any of the other heavy elements that are 
prominent today. 

Indeed, to date not a single star has been observed 
that does not contain metals. Instead stellar nucleosyn- 
thetic products are everywhere. The lowest metallicity 
galaxies known are enriched to ~ 0.02Z© (Searle & Sar- 
gent 1972); the lowest-density regions of the intergalactic 
medium (IGM) appear to be enriched out to the highest 
rcdshifts probed (Schaye et al. 2003; Pettini et al. 2003; 
Aracil et al. 2004); and extensive searches for metal-poor 
halo stars have failed to uncover any PopIII candidates 
beyond a handful of stars with extremely low abundances 
of elements heavier than magnesium (Christlieb et al. 
2002; Cayrel et al. 2004). 

Interpretation of these stellar abundances is especially 
complicated by the fact that searches for metal-poor stars 
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are targeted at the Galactic halo, where dust extinc- 
tion is minimal and crowding is not a serious concern. 
These analyses have provided a number intriguing con- 
straints on the enrichment history of the halo (Freeman 
& Bland-Hawthorn 2002; Beers & Christlieb 2005), in- 
cluding detailed measurements of the metallicity distri- 
bution function from [Fe/H] ~ —5 to [Fe/H] ~ —1 (e.g. 
Ryan & Norris 1991; Beers et al. 1992; Barklem et al. 
2005), a detection of a shift in iron peak ratios below 

[Fe/H] 3 (McWilliam et al. 1995; Carretta et al. 

2002; Cayrel et al. 2004, Cohen et al. 2004), and the 
presence of extremely heavy element deficient stars in 
which [C/Fe] > 2.0 (Christlieb et al. 2002; Frebel et al. 

2005) . While these and other observations have provided 
the basis for several theoretical analyses of the first stars 
(e.g. Hernandez & Ferrara 2001; Scannapieco, & Broad- 
hurst 2001; Oey 2002; Schneider et al. 2003; Tumlinson 

2006) , these are subject to substantial uncertainties due 
to the unknown spatial distribution of such objects. 

In fact, it is still unclear how the observed population 
of halo stars is related to PopIII star formation. Indeed, 
White & Springel (2000) used high-resolution cosmologi- 
cal simulations to show that the oldest stars are naturally 
strongly concentrated towards the center of the Galaxy, 
but stars that form in the full population of T v ; r « 10 4 
K halos are widely distributed. More recently, Diemand, 
Madau, & Moore (2005) used a resimulation technique 
to track the positions of particles contained in the pro- 
togalactic halos collapsing from very rare 3. Oct and 3.5ct 
perturbations (where ct 2 is the variance of linear fluctu- 
ations on the scale of the objects). Associating PopIII 
stars with 3. Oct perturbations, they found that the den- 
sity of such objects in the Galactic halo at the radius of 
the solar orbit is three orders of magnitude lower than 
in the bulge. Associating PopIII stars with 3.5ct per- 
turbations led to even more extreme results, decreasing 
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the number of first stars at the solar orbital radius by 
almost another order of magnitude. Similarly, Karlsson 
(2006) was able to reproduce a "metallicity desert" be- 
tween [Fe/H] = -5 and -4 by associating metal-free star 
formation with only the very earliest forming stars. 

Yet this association of metal-free stars with the high- 
est density regions is likely to be misleading. As shown 
in Scannapieco, Schneider, & Ferrara (2003, hereafter 
SSF03; see also Schneider et al. 2006) cosmological en- 
richment was a local process, which occurred in different 
regions over an extended redshift range. Thus while the 
peak of PopIII star formation is likely to have occurred 
at z ~ 10, such stars should have continued to form 
at much lower redshifts. Furthermore, in these models 
metal free stars do not form in perturbations of a given 
a, but rather within 10 7 - 5 M Q - 10 8 M Q objects, which 
are just large enough to cool within a Hubble time, but 
small enough that they are not clustered near areas of 
previous star formation. 

This spread in a and formation redshift raises the 
prospect that the remnants of metal-free protogalactic 
clouds would have indeed ended up in the halo. This 
would mean that the absence of observed metal-free ob- 
jects might require the first stars to form with an ini- 
tial mass function (IMF) biased to stars with lifetimes 
shorter than a Hubble time, and perhaps even biased to 
masses > IOOMq as suggested by a number of recent the- 
oretical studies (e.g. Nakamura & Umemura 1998; Abel, 
Bryan, & Norman 2000; Bromm et al. 2001; Schneider 
et al. 2002; Ripamonti et al. 2002; Tan & McKee 2004). 
Furthermore, it would mean that some of the stars with 
unusual yields found in the halo of our galaxy are indeed 
providing us with direct information about the yields of 
PopIII stars. 

In this paper we combine an extremely high-resolution 
N-body simulation of the formation of the Milky- Way 
with a semi-analytical model of metal enrichment that 
closely parallels that discussed in SSF03. While we plan 
to conduct a careful comparison of this approach with 
a more direct chemodynamical model of the Milky- Way 
in a companion publication (Brook et al. in prep.), the 
dark-matter based model described here nevertheless has 
a number of distinct advantages. Although gas cooling 
is essential to reproduce the spatial distribution of late- 
forming metal-rich Galactic disk stars, the spatial dis- 
tribution of low-metallicity stars in the bulge and halo 
can be approximated from an N-body point of view. 
Such stars were formed at early times in star clusters 
contained in much smaller dark matter halos, and they 
subsequently evolved as collisionless particles. Thus by 
carefully tagging N-body particles directly in progenitor 
halos as star-particles, we can use them to infer the po- 
sitions of the earliest stellar populations. This approach 
allows us to work at a mass resolution (7.8 x 10 5 M Q per 
dark-matter particle) that is currently impossible in gas 
simulations of the Milky- Way, which means that we can 
capture the formation of all halos within which efficient 
star- formation is expected to occur. 

Our semi- analytic model of metal enrichment, on the 
other hand, gives us the flexibility to use this simulation 
to study a wide range of models for the properties of the 
first stars. As the typical mass of such stars remains un- 
certain within at least two orders of magnitude, the metal 
and kinetic energy output from SNe arising from PopIII 



stars remains completely unknown. This uncertainty is 
further compounded by the possibility that metal-free 
stars in the « 200M© mass range may have given rise to 
tremendously powerful pair-production supernovae (e.g. 
Heger & Woosley 2002), with kinetic energies up to 10 53 
ergs per event. While exploring this full range of possi- 
bilities by direct simulations is extremely computation- 
ally expensive, our technique allows us to carry out this 
parameter study in a few minutes, leading to a number 
of robust conclusions despite the large uncertainties in- 
volved. 

The structure of this work is as follows. In §2 we de- 
scribe our simulations and semi-analytical model. In §3 
we apply our model to consider the distribution of the 
first and second stars in a wide range of cosmological 
scenarios. Conclusions are given in §4. 

2. METHOD 

To identify the location of PopIII stars in our galaxy, 
we combine a high-resolution N-body simulation of the 
formation of a Milky- Way size galaxy with a semi- 
analytical model of metal enrichment. In §2.1 we de- 
scribe our simulation in detail, as well as the method we 
use to extract gravitationally bound halos. In section 
§2.2, we describe our semi-analytical model and discuss 
how it is used to construct the expected positions for 
stars forming from gas that is initially primordial, which 
we label as "the first stars," and stars forming from gas 
that is initially enriched by only primordial stars, which 
we label as "the second stars." 

2.1. Simulation and Identification of Halos 

Throughout our study we assume a Cold Dark Matter 
model with parameters h = 0.71, iloh 2 = 0.135, J7a — 
1 - O , Vbh 2 = 0.0224, cr 8 = 0.9, and n = 1, where h 
is the Hubble constant in units of 100 km/s/Mpc, £! , 
f^A, and Of, are the total matter, vacuum, and baryonic 
densities in units of the critical density, cr| is the variance 
of linear fluctuations on the 8/i _1 Mpc scale, and n is the 
"tilt" of the primordial power spectrum (e.g. Spergel et 
al. 2003). Our N-body simulation is carried out with 
the GCD+ code (Kawata & Gibson 2003a), and it uses a 
multi-resolution technique (Kawata & Gibson 2003b) to 
achieve high-resolution in the regions of interest, while 
the outer regions exerting tidal forces are handled with 
lower resolution. 

The initial conditions for these simulations at z = 
56 are constructed using the public software GRAFIC2 
(Bcrtschingcr 2001). The low-resolution region is 
a 20h~ 1 Mpc diameter sphere selected from a low- 
resolution cosmological simulation, for which an iso- 
lated boundary condition is applied. The highest res- 
olution region is a sphere with a radius four times the 
virial radius of the system at at z = 0. In this re- 
gion, the Dark Matter (DM) particle mass and soften- 
ing length of are M vir = 7.8 x 10 5 M Q and r vir = 540 
pc, respectively. The system consists about 10 6 parti- 
cles within the virial radius. The virial mass and radius 
of the system at z = are M v i r = 7.7 x 10 11 Mq and 
r v ; r = 239 kpc, which are roughly consistent with the 
estimated mass (M v j r = 1.0 x 10 12 M ) and virial ra- 
dius (r vir = 258 kpc) of the Milky Way (Klypin, Zhao, & 
Somerville 2002; Battaglia et al. 2005, see also Dehnen, 
McLaughlin, & Sachania 2006 for a slightly more mas- 
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sive, M v i r = 1.5 x 10 12 M , and compact, r vlr = 200 kpc, 
model.) The initial conditions are essentially the same 
as the KGCD model in Bailin ct al. (2005), but the cur- 
rent simulation has a much higher resolution and slightly 
different cosmological parameters. 

Using a low-resolution simulation including gas physics 
and star formation, we confirm that this initial condition 
will lead to a disk galaxy. A closer comparison between 
our results here and a lower-resolution gas-dynamical 
model will be presented in a future publication (Brook 
et al. in prep.). While gas physics is essential to repro- 
duce the spatial distribution of disk stars, an N-body ap- 
proach is appropriate for the halo and bulge populations 
on which we are focused here. Such stars should initially 
form in much smaller gas disks condensing in neighboring 
halos, and later interact as collisionlcss particles. Thus 
by carefully tagging N-body particles directly in progen- 
itor halos as star-particles, we can use them to infer the 
positions of the first and second stars in our final galaxy. 

The simulation data is output every 0.11 Gyr, and 
at each output, we use a friend-of- friends group finder 
(FOF, Davis et al. 1985) to identify the virialized DM 
halos. We apply a standard approach with a linking pa- 
rameter of b = 0.2 and threshold number of particles of 
50, and we explore the results of varying this parame- 
ter below. We then construct a merger tree, associat- 
ing the halos at different time steps. In particular, we 
search for the halo at each time step n that contains the 
largest number of the member particles of a given halo 
at the previous time step, n — 1. We define the virial 
mass and radius for each halo, taking into account the 
cosmology and redshift, following the fitting formula in 
the Appendix of Kitayama & Suto (1996). 

A key issue is the minimum mass of the halos in which 
stars formed. Before reionization, this depends sensi- 
tively on the presence of H2 {e.g. O'Shea et al. 2005). In 
small objects, molecular hydrogen is easily photodissoci- 
ated by 11.2-13.6 eV photons, to which the universe is 
otherwise transparent. Thus emission from the very first 
stars quickly destroyed all avenues for i? 2 cooling (Dckcl 
& Rees 1987; Haiman, Rees, & Loeb 1997; Ciardi, Fer- 
rara, & Abel 2000). This raised the minimum virial tem- 
perature necessary to cool effectively to approximately 
10 4 K, although the precise value of this transition is the 
subject of debate {e.g. Glover & Brand 2001; Yoshida et 
al. 2003) and is somewhat dependent on the level of the 
high-redshift X-ray background (Haiman, Rees, & Loeb 
1996; Oh 2001; Machacek, Bryan, & Abel 2003). 

A natural question is whether the few stars that dis- 
sociated primordial molecular hydrogen may have been 
also enriched a significant fraction of the intergalactic 
medium. To estimate this fraction we consider two cases, 
one in which the primordial stars that contributed most 
to the dissociating background were « 1OM , and a sec- 
ond case in which the primordial IMF was peaked to form 
very massive stars « 2OOM . In the 10M© case the rate 
of dissociating photons is « 10 47 s _1 (Glover & Brand 
2001; see also Schaercr 2002; Venkatesan, Tumlinson, & 
Shull 2003) and the stars live 30 Myrs (Fagotto et al. 
1994) to generate a total of « 10 62 dissociating photons 
per star. In the 200 M case the dissociating rate is 
w 2 x 10 49 s _1 , but the lifetime is only 2 Myrs, resulting 
in a total of w 10 63 such photons per star. Combining 
these values with the 3 x 10~ 4 mean cosmological frac- 



tion of H2 {e.g. Yoshida et al. 2003), we can estimate 
the typical radius of a region dissociated by a primordial 
star as rdi ss = 200A~ 1//3 comoving kpc in the 10 M case 
and rdi ss = 450A~ 1 / 3 comoving kpc in the 200 M case, 
where A = p/p is the mean overdensity of the dissociated 
region. 

On the other hand, adopting a simple Sedov- Taylor 
estimate of the maximum distance that a given super- 
nova will be able to distribute metals as a function of 
time gives r enrich = 25(1 + z) 2 / 5 ^^ 1 ^- 1 / 5 comov- 
ing kpc where tcyr is the expansion time in Gyrs and 
E 5 i is the energy input in units of 10 51 ergs. If we 
assume no kinetic energy losses, replace icyr with the 
Hubble time, and take a typical redshift of 15, this gives 
Enrich = 35El{ 5 A^ 1 / 5 comoving kpc. 

Since the formation of these very first stars is self- 
truncated by H 2 dissociation, they should cease to form 
and distribute metals when the filling factor of their dis- 
sociation regions reaches unity. This means that the rel- 
ative ratio of the dissociation/metal enrichment volume 
should be an accurate predictor of the total filling fac- 
tor of such very early metals. In the 10 M case this 
ratio is (r onr i c h/fdiss) 3 ~ 0.005 while in the 2OOM case 
(''cnrich/T'diss) 3 ~ 0.002 where we have taken a typical SN 
energy of 10 52 ergs (Heger & Woosley 2002). Thus less 
than 1% of the universe is likely to have been enriched 
by this earliest generation of stars, and this only in the 
very densest regions. Rather the majority of primordial 
enrichment is should have occurred in larger objects with 
T vir > 10 4 K. 

In this case, (sometimes referred to Population II. 5; 
e.g. Grief & Bromm 2006) efficient atomic-line cooling 
establishes a dense locally-stable disk, within which non- 
equilibrium free electrons catalyze i?2- Unlike in less 
massive halos, H2 formation in these objects is largely 
impervious to feedback from external UV fields, due to 
the high densities achieved by atomic cooling (Oh & 
Haiman 2002). As discussed in SSF03, regions of late pri- 
mordial star formation may have already been detected 
as a subgroup of high-redshift Lyman-alpha emitters, al- 
though distinguishing them from Population II/I stars is 
extremely difficult (Dawson et al. 2004). 

In this study then, we adopt a fixed threshold of T vlr > 
10 K for star formation in halos at all redshifts. In our 
cosmology, and for the high redshifts relevant for the 
formation of primordial objects, this corresponds to a 
minimum mass of 

M min = 3.0x 10 9 (l + z)- 3 / 2 M Q . (1) 

While at high-redshift this mass can be over four or- 
ders of magnitude smaller than the final virial mass of 
the Galaxy, our high-resolution simulation nevertheless 
makes it possible to trace these objects. This means that 
at redshifts below 17, all the halos relevant to our study 
will contain at least 50 particles and be well-identified by 
our FOF group finder. However, at the very earliest red- 
shifts > 17 such halos will only be marginally resolved. 

2.2. Identification of First and Second Objects 

Having constructed the history of dark matter halos 
as a function of redshift, our next step is to use this to 
identify two types of objects: i) halos that collapse out 
of primordial gas, which we identify as Population III 
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objects containing "the first stars," and ii) halos that 
collapse from gas that has been enriched purely by ma- 
terial ejected from the Population III objects, which we 
identify as second-generation objects containing "the sec- 
ond stars." While we vary the IMF of the first stars over 
a wide range of possibilities, we assume that the sec- 
ond stars form with an IMF similar to that seen today, 
although they are likely to display peculiar enrichment 
patterns. 

It is important to make clear that our model is not 
able to say anything about the progression of star forma- 
tion and self-enrichment within a particular halo (see e.g. 
Susa & Umemura 2006). Rather, our labels purely reflect 
the properties of the most metal-poor stars formed in any 
given region. Thus it is entirely possible that some frac- 
tion of the stars forming in halos that collapse from pri- 
mordial gas will nevertheless contain metals that were in- 
herited from other stars forming in the same temporally- 
extended burst. For our purposes here then, stars formed 
in such self-enriched primordial star clusters are included 
in the distribution of "the first stars." Similarly we ex- 
clude all halos that are direct descendents of first objects 
from our list of second generation objects, as these would 
exist cospatially with stars that were truly primordial. 

These definitions are intended to extract the most in- 
formation possible from our modeling approach. For ob- 
servers, anywhere labeled first stars would be a good 
place to search for metal-free stars, if they are sufficiency 
long-lived to survive until today. Similarly, both regions 
labeled first stars and second stars in out method would 
be a good places to search for stars enriched purely by 
the products of metal-free stars. From a theoretical point 
of view, the spatial distribution of first stars in our mod- 
els should be directly comparable to the spatial distri- 
bution of metal-free stars in future detailed gas simula- 
tion, while the spatial distribution of both first and second 
stars should be comparable to the spatial distribution of 
PopIII-enriched stars in such simulations. These defini- 
tions are important to keep in mind when interpreting 
our results, as emphasized below. 

As the distribution of both these populations of stars is 
primarily dependent on the highly-uncertain properties 
of Population III objects, our goal is to develop a simple 
model of metal dispersal that allows us explore a large 
range of possibilities in a straightforward way. Following 
our approach in SSF03, we adopt a model of outflows as 
spherical shells expanding into the Hubble flow (Ostriker 
& McKee 1988; Tegmark, Silk & Evrard 1993), for both 
PopIII and PopII/I objects. These shells are assumed 
to be driven only by their internal hot gas pressure and 
decelerated by inertia due to accreting material and grav- 
itational drag while escaping from the host. In this case 
there are two relevant evolutionary equations: one that 
describes the change in the velocity of the bubble, 
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where the terms on the right hand side correspond to 
internal pressure, accumulation of material from the sur- 
rounding medium, and gravitational drag, respectively, 
and a second equation that describes the evolution of 
the internal energy of the bubble, 



E b = L(t) - 4irR 2 s R s P b , 



(3) 



where the terms on the right hand side correspond to 
energy input from the supernovae and the PdV work 
done by the interior gas. In these equations the over- 
dots represent time derivatives, the subscripts s and b 
indicate shell and bubble quantities respectively, R s is 
physical radius of the shell, E b is the internal energy of 
the hot bubble gas, P b is the pressure of this gas, and p 
is the mean IGM background density. Here we assume 
adiabatic expansion with an index 7 = 5/3 such that 
P b = E b /2nR 3 s . 

The evolution of each such bubble is completely de- 
termined by the mechanical luminosity evolution as- 
signed, L(t). Again for simplicity, we take both PopIII 
and PopII/I objects to undergo a starburst phase of 
isN = 5 x 10 7 years, but with different prescriptions for 
the total energy input into the wind. We approximately 
account for the gravitational potential of the host galaxy 
by subtracting the value of GM 2 (Vl b /Vtm)/r Y - a from the 
total wind energy, where M is the total mass (dark + 
baryonic) of the object and r v ; r is its virial radius. Ap- 
proximating the collapse overdensity by a constant value 
of 178, the overdensity associated with a virialized cos- 
mological object, then gives a mechanical luminosity of 

L(t) = 16OL [fl^f.E^N 111 ' 11 (4) 

-5 x 10- 12 (A/ fc /ift M M) 2/3 (l + *c)] M b e{t SN - t), 

where /* is the fraction of gas converted into stars, f w 
is the fraction of the SN kinetic energy that is channeled 
into the galaxy outflow, E^ n in the kinetic energy per 
supernova, Af is the number of supernovae per M Q of 
stars, M b is the baryonic mass of the galaxy in units of 
solar mass, z c is the collapse redshift of the object, and 
9 is the Heaviside step function. Note the precise value 
of isN is not important as long as it short relative to the 
« 500 Myr Hubble time at the redshift in which we are 
most interested. However, the value of the prefactor in 
front of £"kin in eq. ^ is defined by this choice. Halving 
isN, for example would require for us to double this factor 
but have only a very weak effect on the time that the 
bubbles reach their neighbors, and hence on the final 
distribution of first and second stars. Note also that the 
gravitational drag is only important in objects with total 
masses > 10 11 (1 + z)~ s / 2 Mq, and thus has only a minor 
impact on our results here. 

The only difference in the wind evolution of 
PopIII and PopII/I arises from the product, 
fl U,U f w E^' u N 111 ' 11 , which, again following SSF03, 
we define as the "energy input per unit gas mass" 
Sg 11 ' 11 , which we express below in units of 10 51 ergM Q " 1 . 
To determine this value in the PopII/I case, we take f* 1 
to be 0.1, which gives good agreement with the observed 
high-redshift star formation rates and abundances of 
metals measured in quasar absorption line systems 
(Thacker, Scannapieco, & Davis 2002; Scannapieco, 
Ferrara, & Madau 2002). Also in this case, we constrain 
f w by combining the overall efficiency of 30% derived 
for the 2 x 10 8 M Q object simulated by Mori, Ferrara, & 
Madau (2002) with the mass scaling derived in Ferrara, 
Pcttini & Shchckinov (2000), which was obtained by 
determining the fraction of starburst sites that can 
produce a blow-out in a galaxy of a given mass. This 
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gives f w (M) = 0.38 B (M)/5 B (M = 2 x 10 8 M Q ) where 

{ 1.0 N t < 1 

S B (M) = I 1.0 - 0.165 ln(N t ) 1 < A> t < 100 , 

[ [1.0 - 0.165 ln(100)] 100 jV t -1 100 < N t 

( 5 ) 

and Nt = 1.7 x 10 7 ' (SIi,/Q,m)M /Mq is a dimensionless 
parameter that increases linearly with the overall num- 
ber of SNe produced in a starburst, divided by the star 
formation efficiency, f* . 

In the PopIII case, on the other hand, there are no di- 
rect constraints on either f£ n or the wind efficiency. For 
these objects we allow these parameters to be free, vary- 
ing Eg 11 over a large range as discussed below. Finally, 
when outflows slow down to the point that they are no 
longer supersonic, our approximations break down, and 
the shell is likely to be fragmented by random motions. 
At this point we let the bubble expand with the Hubble 
flow. These prescriptions then give us a complete model 
of Rhiji(M, z, z'), the comoving radius of an outflow at 
a redshift z' associated with a PopIII or PopII/I halo of 
mass M collapsing at a redshift z. 

Our goal is then to tag halos that form outside of all 
expanding shells as first stars objects, and to tag halos 
that form within PopIII outflows as second-stars objects. 
In our fiducial model we assume that in each such halo 
stars form with the same radial profile as the dark mat- 
ter, although we also vary this assumption in §2.2. As 
we are concerned only with whether a halo is contained 
within any such wind from a neighbor or direct progen- 
itor, we need not worry about the possibility of double 
counting by identifying the same collapsed object at mul- 
tiple outputs. Rather to determine if a halo forming at 
a redshift z' has been enriched, we can safely compute 
Riii,ii(M, z, z') for each and every object forming with 
z > z' and determine if any of them has overtaken the 
lower- redshift object. In this case the impact of a given 
halo on its neighbors as well as its direct descendants 
will be determined by the maximum Riuji[M(z), z, z'\ 
computed at any z > z' , while all other smaller shells 
will have no impact on the calculation. 

While this picture neglects the precise density and in- 
fall structure around each forming object, we are able 
to quantify the errors arising from this approximation 
by adopting two methods for computing the distance be- 
tween two objects. Note however that as the outflows we 
are most concerned with occur at very high redshift, the 
nonlinear evolution of structure will not be as severe an 
issue as it is for lower-redshift winds. 

One method of computing distances, which we shall re- 
fer to as "Eulerian," allows each wind to move along with 
the center of mass of its central halo. This means that as 
we move between each output we compute the most likely 
daughter halo associated with each higher redshift out- 
flow and re-center it around the position of this object. 
This takes advantage of the nonlinear motions computed 
in the simulation, but is likely to overestimate the impact 
of winds as it does not take into account the impact of 
infall around collapsing structures. 

A second method of computing distances, which we 
shall refer to as "Lagrangian," works completely in the 
reference frame in which the position of each halo is com- 
puted as the center of the mass of the original z = 55.7 
positions of the particles it contains. While this ap- 
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proach neglects much of the nonlinear information avail- 
able from the simulation, it is nevertheless likely to rep- 
resent a more fair comparison as: i.) our thin-shell so- 
lution, which assumes expansion in the Hubble flow, is 
essentially already working in this coordinate system; ii.) 
enrichment should be more closely dependent on the col- 
umn depth between the source and recipient halos than 
on their physical separation. Furthermore this approach 
most directly parallels that of our previous analytical 
models (SSF03), and thus we adopt it as our fiducial 
method in the results presented below. 

3. RESULTS 

In this section we summarize our results for a wide 
range of parameters choices and assumptions, the details 
of which is described in detail below. For reference pur- 
poses, the specifics of each of these runs is summarized 
in Table 1. 

3.1. Fiducial Modeling 

In Figure [TJ we show the star formation rate in first 
and second stars as a function of redshift, for a wide 
range of feedback efficiencies and for both Lagrangian 
and Eulerian distances, corresponding to runs 2L-4L and 
2E-4E in Table 1 respectively. When computing the star 
formation rate, we consider only stars that end up grav- 
itationally bound to our final galaxy, so as to exclude 
objects formed near the boundary of the simulation vol- 
ume. In all cases we assume fl 11 = fj 1 = 0.1 such that 
10% of the gas mass in each PopIII halo is converted 
into stars. Note that while this value is consistent with 
the observed star formation rate of normal stars at lower 
rcdshifts (e.g. Scannapieco, Ferrara, & Madau 02), it is 
certainly possible that metal-free star formation was less 
efficient and that a lower value such as // 7/ = 0.01 is 
appropriate in the PopIII case. While this normalization 
will remain a primary uncertainty in our investigation, it 
can be easily accounted for by rescaling our final results, 
and does not affect the overall trends on which we are 
primarily focused. 

In SSF03, we considered the range of physically pos- 
sible Eg 11 values in detail. An upper bound is pro- 
vided by the case in which the initial mass function 
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Fig. 1. — Star formation rate as a function of redshift in our 
models. In each panel, the solid lines show the rate in first stars 
objects containing initially primordial gas, while the dashed lines 
show the star-formation rate in second-stars objects, which have 
been pre-enriched by only the products of primordial stars. From 
left to right the panels represent models ranging from those con- 
taining extremely strong (£g = 10 -2 ) winds, to those with much 
more gradual (Si 1 = 10 — 4 ) metal ejection. A 10% star formation 
efficiency is assumed in all cases. The top row shows the results of 
models in which the distances between halos have been computed 
using the "Lagrangian method" described in §2.2, and the bottom 
row represents models in which distances have been computed us- 
ing the "Eulerian method," which predicts a somewhat stronger 
impact from winds. Each panel is labeled by its corresponding 
model name as described in Table 1. Despite the large range of 
parameters and methods considered, metal-free star formation oc- 
curs down to z 3 in many cases, and down to z ss 5 even in the 
most extreme estimates. 

(IMF) of primordial stars is biased to very high masses 
> IOOMq. While much of the stars in this mass range 
would have died forming black holes, progenitors with 
masses between 140 and 260M Q end their lives in tremen- 
dously powerful pair-production supernovae. In these 
SNe, e + e~ pair creation softens the equation state at the 
end of central carbon burning, leading to rapid collapse, 
followed by explosive burning (e.g. Barkat, Rakavy, & 
Sack 1967; Ober, El Eid, & Fricke 1983; Bond, Arnett, 
Carr 1984; Heger & Woosley 2002). The more massive 
the star, the higher the temperature at bounce and the 
heavier the elements that are produced by nuclear fusion. 
In all cases, the star is completely disrupted and the ejec- 
tion energies (« 3 — 100 x 10 51 ergs) are enormous. In this 
case, choosing a narrow Gaussian IMF centered at 200 
Mq, taking a 10% star formation efficiency, and fixing 
f w « 0.5 places a generous upper limit of E^ 11 = 10~ 2 , 
corresponding to the strongest feedback model in Fig. [TJ 
Note that while this extreme model is useful to set an up- 
per limit on feedback, widespread enrichment from stars 
with such an IMF is inconsistent with the observed abun- 
dance patterns of metal poor stars owing to the absence 
of pair-instability yields in the data (Umeda & Nomoto 
2003; Tumlinson, Venktesan, & Shull 2004). 

A lower limit on Eg 11 , on the other hand, is computed 
by assuming that such stars form with a typical Salpeter 
IMF, fl u « 0.1 and f w « 0.1, which gives E^ H « 10~ 4 . 
Not only does this model take metal ejection from PopIII 
to be somewhat less efficient than from PopII/I objects, 



but as shown in SSF03, even less efficient models of 
PopIII metal ejection are so weak that a large primor- 
dial star formation rate is predicted to this day. Thus 
the weakest, Eg 11 = 10~ 4 , model shown in Figure [T] can 
be taken to represent a conservative estimate of the min- 
imum efficiency with which metals from PopIII objects 
were mixed into their surroundings. 

This large range of efficiencies of PopIII objects is what 
forces us to adopt a semi-analytical model of winds in 
the first place. Yet, as shown in Figure [IJ robust gen- 
eral conclusions can be drawn despite these uncertainties. 
Firstly, the primordial objects that predated the Milky- 
Way were formed over an extremely large range of red- 
shifts, mirroring the overall cosmological PopIII distribu- 
tion (sec SSF03). Furthermore, for all Eg 11 values, and 
in both the Eulerian and Lagrangian models, appreciable 
PopIII star formation continues down to a relatively low 
redshift of z ~ 5. In fact, excluding the most extreme 
models, PopIII formation continues at a substantial rate 
down to z ss 3. 

While increasing the efficiency of feedback has the im- 
pact of moving the PopIII star formation to higher red- 
shifts, the overall z » 10 peak and Az w 5 width of 
the epoch of PopIII formation are surprisingly constant 
among models. Instead, the distribution of second stars 
is much more sensitive to these uncertainties, ranging 
from an extended burst in the Eg 11 = 10~ 2 case, to a 

weak high-redshift burst in the Eg 11 = 10~ 4 case. 

Again, we stress that these second stars are the stars 
forming in halos that have been initally pre-enriched by 
PopIII neighbors, and that a significant population of 
Pop-Ill enriched stars are also likely to exist in "first 
stars" halos forming of initially primordial gas. Finally, 
comparing the Eulerian and Lagrangian models, we find, 
as expected, that Eulerian distances lead to somewhat 
more efficient metal dispersal, although this is likely to 
be somewhat of an overestimate. Nevertheless the differ- 
ences between these two models are relatively mild for 
all choices of Eg 11 , giving us further confidence in the 
approximations involved in our semi-analytic approach. 

Figure[5]demonstrates the positions of the first and sec- 
ond stars at different redshifts in the Lagrangian model 
with E^ 11 = 10~ 2 . Here and below we compute the 
masses as they would be if there were no mass loss due 
to stellar evolution. Thus these plots should be rescaled 
by the fraction of first and second stars surviving to a 
given redshift, model-dependent quantities that we do 
not address here. 

At early times (z — 9.84), the first stars are forming 
close to the central density peak of the progenitor galaxy, 
due to the higher density peaks in this region. Second 
stars form in the halo in the neighborhood of first stars, 
because they condense from gas that is enriched by the 
material from the explosions from these objects. At a 
later time (z = 6.02), new first stars are still forming, 
but now on the outer regions of the progenitor galaxy, 
because they are not yet affected by the bubbles from 
the first stars that formed preferentially in the central 
region. As seen in the star formation rate plot, Figure 
[U the formation of the first and second stars is complete 
around z = 3, at which time the full region is enriched 
with metals. Also at this redshift, first and second stars 
start to be accreted into the assembling Galactic halo. 
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Fig. 2. — Distributions of first and second stars (3rd and 4th columns) in a a 200 proper kpc a region, at various rcdshifts, in model 2L, 
a Lagrangian model with £g U = 10~ 2 . For comparison, the 1st and 2nd columns show all the particles and all the particles within halos 
with masses above M m j n in equation (1), respectively. In the 1st and 2nd columns, to avoid confusion, we only plot 1/10 of the whole 
sample. 
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Fig. 3. — Distributions of first and second stars at z = in our fiducial, Lagrangian models. From top to bottom, the rows correspond to 
the central 100 kpc region in the £g H = 10~ 2 model, the central 100 kpc in the £g 1 = 10~ 4 model, the central 10 kpc in the S^ 11 = 10~ 2 
model, and the central 10 kpc in the £^ n = 10~ 4 model. In interpreting this figure, it is important to note that we are explicitly ignoring 
mass loss due to stellar evolution. 
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Fig. 4. — Radial profile of first stars forming from initially pri- 
mordial gas (left column) and second stars formed from initially 
PopIII enriched gas (right column) at z = in our fiducial La- 
grangian models. In the top row, the dotted lines give the overall 
dark matter profile of the galaxy, which is compared with the radial 
density of first and second stars in models 2L (solid), 3L (short- 
dashed), and 4L (dot-dashed). The lower row shows the fraction 
of the total density in first and second stars, with symbols as in 
the upper panels. Like Fig. \3\ and those below, this plot does not 
include any mass loss due to stellar evolution. Furthermore, it 
assumes f* 1 1 = fl 1 = 0.1. 



This assembly is almost complete by z = 1, and thus the 
distribution at this redshift is similar to that at z = 0, al- 
though at z = the stars are more smoothly distributed. 

In Figure [3] we show the z = spatial distribution 
of the first and second stars in closer detail, contrast- 
ing two extreme models with strong £^ n = 1CP 2 and 

weak Sg 11 = 10~ 4 feedback. From this figure, we see 
that for all choices of S^ 11 , the first or second stars are 
not confined to the center of the Galaxy. Rather, this 
plots shows a widespread distribution that is punctuated 
by dense concentrations of stars associated with satel- 
lite galaxies, as seen most clearly in the 100 kpc pan- 
els. Furthermore, comparing the satellite distributions 
between the two models shows that as one moves from 
the strongest to the weakest feedback models, a large 
number of the satellite second-generation objects switch 
over to the first-generation population. This is because 
these objects were only enriched by primordial stars, thus 
a large number of them form from primordial gas if the 
winds from PopIII stars are weakened. The more iso- 
lated the satellite galaxy, the more likely it is to contain 
primordial stars even in strong feedback models. 

Focusing on the central 10 kpc of our simulation, which 
is more representative of the Galaxy itself, we see that 
the distribution of first and second stars, while peaked 
towards the center, nevertheless contains many stars at 
large radii. Again this is true even in the most extreme 
feedback models, and it represents the single most impor- 
tant finding of our investigation. If they have sufficiently 
long lifetimes, a significant number of stars formed in 
initially primordial star clusters should be found in the 
Galactic halo, regardless of the specifics of PopIII star 
formation. 

To quantify this statement further, in Figure|4]we plot 



the radial mass density of first and second stars, as com- 
pared to the dark matter distribution in our simulation. 
Again, in computing these masses, we assume that 10% 
of the gas mass in each halo is converted into stars and 
we intentionally make no attempt to account for mass 
loss due to the evolution of these stars. Furthermore, we 
plot only out to 250 kpc, which is 1/4 the size of the 
high-resolution region. 

Here we see that the density profiles of first and second 
stars are similar to the total dark matter density profile, 
although the second stars have a slightly shallower slope. 
As a result, the density of first stars at the center is 100 
times higher than at the 8 kpc orbital radius of the sun. 
However, the important number for developing observa- 
tional strategies is the relative density of such stars with 
respect to field stars. While this is not directly computed 
in our simulation, the lower panels of Figure 2] show the 
local mass density of stars normalized by the local den- 
sity of dark matter. 

Amazingly, the mass fraction contained in the first 
stars varies only very weakly with radius. Moving from 
1 to 100 kpc in the S^ 11 = 10~ 4 model, for example, 
the fraction of the mass in first stars decreases only by 
a factor m 4. Furthermore, increasing the efficiency of 
PopIII winds to £^ n = 10 -2 has the effect of decreasing 
the fraction of first stars without strongly affecting their 
radial distribution. Note that in the inner region, the 
overall stellar density is expected to be higher than dark 
matter density (e.g. Widrow & Dubinski 2005), that is 
the density of spheroidal stellar component of the Milky 
Way decreases with the radius more quickly than the 
dark matter component. Therefore, the fraction of the 
first stars compared with the field stars should be some- 
what higher at large radii than what Figure |4] predicts. 
This raises the possibility that the overall fraction of first 
stars may even be higher in the halo than in the bulge! 

These results have a number of immediate implica- 
tions. First, as stars enriched only by primordial stars 
(which trace the overall distribution of both first and sec- 
ond stars in our study) are found throughout the Galac- 
tic halo, measured abundances of metal-poor stars in the 
Galactic halo should be taken as directly constraining the 
properties of primordial objects. Furthermore since only 
a slightly larger fraction of stars in the Galactic bulge 
are primordial, there is no compelling theoretical reason 
that observations should have to focus in this extremely 
difficult region. Instead the lack of metal-free stars ob- 
served in the halo is likely to imply a real lower mass 
limit in the metal free IMF of at least 0.8M Q , the mass 
of a low-metallicity star with a lifetime comparable to the 
Hubble time (Fagotto et al. 1994). In fact the ~ 1O~ 6 M 
pc -3 value for PopIII stars in the solar neighborhood we 
compute is so high as compared to the observed stellar 
mass density of 5 x 10 _5 Mq pc -3 (Preston, Shectman, & 
Beers 1991) that several primordial 0.8 Afg stars would 
have been observed even if fl 11 were over an order of 
magnitude lower than the 0.1 value we used to normal- 
ize our approach. Finally, it means that there is a good 
chance that many of the stars with extremely unusual 
abundances, such as HE0107-5240 (Christlieb et al. 2002) 
and HE1327-2326 (Frebel et al. 2005; Aoki et al. 2006), 
may in fact represent the products of self-enrichment in 
primordial star clusters, which display the yields of a 
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Fig. 5. — Top: Average formation redshift of first and second 
stars in our Lagrangian models. As in Fig. [4] models L2, L3, 
and L4 are shown by the solid, short-dashed, and dot-dashed lines 
respectively. Bottom: Normalized number density of first stars 
forming above z = 10 and second stars forming above z = 8, with 
lines as in the upper panels. For comparison these plots also include 
the total mass fraction of first and second stars, given by the dotted 
lines. 

single neighboring PopIII star (e.g. Umeda & Nomoto 
2003; Suda et al. 2004; Iwamoto et al. 2005; Tsujimoto 
& Shigeyama 2006). Note that this widespread PopIII 
star distribution is also consistent with studies of the 
lithium-7 abundance of metal-poor [Fe/H] < -1.3 stars, 
which suggest that a large fraction of the baryonic mat- 
ter in the early halo may have been processed through 
such stars (Piau et al. 2006). 

In the upper left panel of Figure [5] we plot the mean 
redshift of formation of the first stars as a function of 
radius. While this is not a directly observable quantity, 
it nevertheless helps in contrasting our results with pre- 
vious estimates that assume a single formation redshift 
or a for the first stars. In all models (zForm) decreases 
monotonically with radius, moving from (zp rm) ~ 13 
near the Galactic center, to (zForm) < 8 at a distance of 
100 kpc. In the lower left panel of this figure, we plot 
the normalized radial mass density of first stars with for- 
mation redshifts above 10. This shows that essentially 
all of the earliest forming first stars are located near the 
Galactic center, and that their relative number density 
drops by at least a factor of 30 as we move towards 8 
kpc. 

These results are reminiscent of those of White & 
Springel (2000) and Diemand, Madau, & Moore (2005), 
who used N-body simulations to show that the oldest 
stars are very strongly concentrated near the Galactic 
center (see also Miralda-Escude 2000). Interestingly, 
White & Springel (2000) also estimated the history of 
very metal-poor Galactic stars by rescaling an N-body 
cluster simulation to Milky Way scales, and associating 
PopIII stars with all particles collapsing into T vu = 10 4 K 
halos. They found that these stars had a mean forma- 
tion redshift of z ~ 4 and a spatial distribution closely 
following the overall dark matter distribution. Our more 
detailed modeling shows that metal free star formation 
displays elements of both these scenarios. Primordial 



Fig. 6. — Radial profile of first and second stars at z = in 
Eulerian models 2E7solid), 3E (dashed), and 4E (dot-dashed). 
Panels are as in Fig [4] Compared to the Lagrangian models, the 
Eulerian models show a somewhat smaller number of first stars and 
a more strongly centrally-peaked concentration of second stars. 

stars form in the smallest halos that can cool efficiently, 
but only within a subset of these halos that is biased 
to the highest redshifts. These stars do not follow the 
overall mass distribution, but display a only weak radial 
density gradient. And their overall formation redshifts, 
while always high, nevertheless drop off strongly as a 
function of radius. 

Finally, the upper right panel of Figure [5] shows the av- 
erage formation redshift of second stars in our simulation. 
Again we remind the reader that these are the stars that 
contain metal from primordial stars but are not formed 
in the same star clusters. These stars display a similar 
radial trend as the first stars, but shifted to somewhat 
lower formation redshifts. In this case, the majority of 
these stars in the solar neighborhood are formed at red- 
shifts below 8, as shown in the lower right panel of this 
figure. 

3.2. Model Tests and Refinements 

As our results in Figs. [3] and 2] have a number of wide- 
reaching implications, it is important to understand to 
what degree they are dependent on our assumptions. 
The primary simplification in our model is our use of 
a Lagrangian coordinate system combined with a bubble 
model that assumes expansion into the Hubble flow. To 
asses this approximation, we show in Figure [6] the results 
of a calculation that replaces the initial coordinates with 
"Eulerian positions" that re-center each bubble from out- 
put to output, as described in §2.2. As this approach 
combines the large bubble radii computed in the absence 
of infall, with the small separation computed by account- 
ing for the proper motions of halos, it can be taken as 
providing a lower limit on the distribution of first stars. 

Thus, it is not surprising that the densities of first stars 
in these models are lower than their Lagrangian counter- 
parts. Yet, in no cases are these differences dramatic. In- 
stead, the radial distributions in the S^ 11 — 10 -3 (model 
E3 in Table 1) and 10 -4 (model E4) cases in Figure 
[4] show only a slightly lower overall normalization and 
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Fig. 7. — Comparison of the z = distribution of first and second 
stars in our Lagrangian and Eulerian models. The top two rows 
compare the results of strong feedback models (£g = 10 — 2 ) with 
a Lagrangian distance (top row) and an Eulerian distance (second 
row). The bottom two rows show a similar comparison for weak 



feedback models (Eg 



10 ) with Lagrangian (third row) and 



Eulerian (bottom row) distances. 

a marginally steeper radial trend than seen in the La- 
grangian case. The £^ n = 10~ 2 model also shows an 
overall radial trend that is very similar to that in Fig. 01 
although it contains a gap from 5-20 kpc with a factor of 
» 2 deficit in stars. Thus, even in this most extreme of 
any of our models, the Galactic halo contains a substan- 
tial number of the remnants of PopIII objects. 

On the other hand, much larger differences between 
our Eulerian and Lagrangian models are visible in the 
distribution of second stars. These stars are substantially 
more centrally-concentrated in the Eulerian models, and 
in the £^ n = 1CP 4 case this difference is so severe that 
very few second stars are located in the Galactic halo. 
This represents the largest uncertainty in our approach, 
and thus for the limited case of a weak-feedback model, 
we are unable to draw any clear conclusions as to the 
halo distribution of second stars. 

The differences between the Lagrangian and Eulerian 
approaches are further investigated in Figure which 
shows the spatial distribution of first and second stars in 
both models. Here we see the same overall trends as in 
the radial profiles. In the strong feedback case, first and 
second stars are more concentrated than in the Eulerian 
approach, although in both models significant numbers of 
these stars are found in the halo and in satellite galaxies. 
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Fig. 8. — Radial profile of first and second stars at z = in our 
Lagrangian models, but now only associating stars with the mass 
in the central 10% of each halo. Panels are as in Figs. [4] and [6] 
while the lines represent models 2Lc (solid), 3Lc (dashed), and 4Lc 
(dot-dashed). In all cases, the radial profiles are consistent with 
those in our fiducial approach, shown in Fig. [4] 

In the weak feedback case, on the other hand, the differ- 
ences in the first-star distribution are minor, but the lack 
of second stars in the Galactic halo in the Eulerian model 
is considerable. Finally, apart from the difference in ra- 
dial trends, the detailed distribution of satellites contain- 
ing PopIII stars differ somewhat between each of these 
approaches. Although in general, there are fewer PopIII 
satellites in the Eulerian cases, there are also occasion- 
ally PopIII satellites in these models without Lagrangian 
counterparts. 

A second important issue is our association of stars 
with the overall distribution of dark matter in each halo. 
In reality, the gas in a given object will condense into a 
rotationally-supported disk whose radius is substantially 
smaller that its dark matter virial radius. Furthermore, 
due to the increased impact of tidal effects on the most 
loosely bound particles, it is possible that stars forming 
in such disks may end up distributed significantly differ- 
ently at z = than their associated dark matter. 

To account for this possibility, we repeated our La- 
grangian model, but now only associating stars with the 
inner 10% of the dark matter particles in first and second- 
generation objects, while leaving the number of stars 
formed in each halo fixed. The resulting radial distribu- 
tions and spatial profiles are shown in Figures [8] and [9l 
respectively. Here we see that while limiting the number 
of particles associated with stars leads to larger statisti- 
cal fluctuations, the distribution of first and second stars 
otherwise remains the same as in our fiducial modeling. 
Carrying out the same exercise in the Eulerian case yields 
similar agreement, and thus our results do not seem to 
be affected by differences in tides across halos. 

Another important parameter in our analysis is our 
choice of b = 0.2 linking length when identifying halos. 
To test the dependence of our results on this choice, we 
generated new lists using a much smaller 0.15 length, 
which associates halos with objects almost twice as dense 
as in the standard picture. In Figure [TU] we show the 
results of this calculation in the Lagrangian case (models 
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Fig. 9. — Comparison of the z = distribution of first and second 
stars between our fiducial models and Lagrangian models in which 
stars are associated with the mass in the inner 10% of each halo. 
Consistent with Fig. [8] there are no strong difference between these 
two approaches. 

2L15-4L15). 

From the left panels of this figure, we see that while 
adopting this more stringent criteria reduces the number 
of PopIII stars somewhat, this effect is no stronger than 
a factor of « 2. Furthermore, the overall radial profiles 
in the b — 0.15 case are very similar, if not flatter than in 
the fiducial model, indicating that the choice of linking 
length has no impact on our conclusions about the radial 
distribution of first stars. On the other hand, changing 
b results in a substantial increase the number of second- 
stars within 8 kpc as shown in the right panels. This can 
be traced to the presence of a large object at z = 11, 
which was identified as a single large halo in the b = 0.2 
case, but split into several smaller objects in the b = 0.15 
case. While in the fiducial case, this large object was a 
direct descendant of a first stars halo, and thus excluded 
from the second-stars count, in the 0.15 case, most of the 
mass in this object was tagged as neighboring halos, and 
was added to the second-stars count. While this differ- 
ence is consistent with our approach to counting first and 
second stars, it nevertheless illustrates some of the am- 
biguities involved in distinguishing stars that are formed 
in gas initially enriched by primordial stars, from stars 
that are formed within purely self-enriched primordial 
gas. Note however, that at Galactocentric radii > 8 kpc, 
the overall distributions of second stars are similar be- 
tween the models with different choices of linking length. 
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Fig. 10. — Radial profile of first and second stars at z = in 
Lagrangian models with a 0.15 choice of linking length. Panels are 
as in figure [4] and [6] and the lines represent models 2L15 (solid), 
3L15 (dashed), and 4L15 (dot-dashed). For the first stars, the 
radial profiles are consistent with those in our fiducial, b = 0.2 
approach, shown in Fig. [4] However, choosing b = 0.15 results 
in more second-stars being identified within « 8 kpc, which is 
primarily due to a large high-redshift object that is split into several 
sub-halos with this smaller choice of linking length. 



Another concern is the possibility that our models are 
assigning first and second stars to all halos that collapse 
in primordial areas or areas enriched only by primordial 
stars. However, a one-to-one association between dark 
matter halos and Milky- Way satellites leads to a large 
excess with respect to their observed numbers, which is 
commonly referred to as the missing satellite problem 
(e.g. Moore et al. 1999). Although ram pressure strip- 
ping before star formation (e.g. Scannapieco, Ferrara, 
& Broadhurst 2002) tidal stripping after star formation 
(e.g. Kravstov, Gnedin, & Klypin 2004) and other issues 
are important in resolving this issue, the largest source 
of suppression of dwarf galaxy formation is the increase 
in IGM thermal pressure following reionization (e.g. . 
Bullock, Kravstov, & Weinberg 2000). While reioniza- 
tion will clearly have the largest effect on the very lat- 
est forming stars, rather that the objects in which we 
are most concerned here, the relatively large number of 
satellites seen in plots such as Figs.[2]and[3]indicates that 
our results may nevertheless suffer from such a excess of 
Galactic satellites. 

Recently Benson et al. (2002), have shown that incor- 
porating the increase in IGM pressure associated with 
reionization into semi-analytical models of galaxy for- 
mation helps to bring their faint end galaxy counts into 
rough agreement with observations. Following their lead 
then, we carried out a set of reionization runs, in which 
we imposed the minimum filtering mass which approx- 
imates that computed in Gnedin (2000). In this case 
the minimum virial temperature in our calculation is in- 
creased linearly from the T vu > 10 4 K molecular hydro- 
gen limit at z = 9, to T v ; r > 5 x 10 4 below the final 
overlap redshift of z = 7. The resulting spatial distribu- 
tions and radial profiles are shown in Figures [Tl] and [T2l 
respectively. 

From these plots, it is clear that incorporating such 
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Fig. 11. — Comparison between the spatial distribution of stars 
in models with a fixed minimum virial temperature (2Lc and 4Lc), 
and models with a gradually increasing threshold associated with 
reionization (2LRc and 4LRc). While incorporating reionization 
leads to a far smaller number of first and second satellites, the 
overall distribution of such stars remains widespread. 
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Fig. 12. — Radial profile of first and second stars at z = 
in Lagrangian models that include reionization. Panels are as in 
Figures [4] and [6] and the lines represent models 2LRc (solid), 3LRc 
(dashed), and 4LRc (dot-dashed). Increasing the filtering scale 
from 10 4 K at z = 9 to 5 X 10 5 K at z = 7 has only a weak effect on 
our results. 
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Fig. 13. — The anisotropy parameter of the velocity dispersion, 
f}, as a function of radius, for all the particles (dotted), first star 
particles (solid), and second star particles (dashed). 

suppression leads to far fewer first and second satellites, 
moving from « 25 such objects with 100 kpc in our fidu- 
cial models, to less than 10 in the reionization models. 
Nevertheless, as shown in Figure[T2]the overall radial pro- 
file of first and second stars are quite similar to our fidu- 
cial results. With metal-free stars (first stars), and stars 
that contain the products of metal-free stars (first + sec- 
ond stars) forming at all Galactocentric radii. Thus we 
do not expect the details of suppression of dwarf galaxy 
formation to have a strong effect on our conclusions. 

Finally, we note that the cosmological parameters used 
in our model have recently been revised by the 3-year 
data release from the Wilkinson Microwave Anisotropy 
Probe experiment (Spergel et al. 2006). In general, we ex- 
pect that the low as and n values implied by these mea- 
surements should strenghten our conclusions. Indeed, 
the reduced power on small scales leads to a delay in 
the onset of star formation and metal enrichment in our 
T v ; r > 10 4 K halos. As a consequence, the first stars will 
form at more recent epochs than in our fiducial model, 
and are likely to give a larger contribution to the halo 
population. 

3.3. Kinematics of First and Second Stars 

Our simulation results offer not only the spatial distri- 
bution of first and second stars, but also their kinematics. 
In this section, we present the results of a representative 
Lagrangian model, model 2Lc. Figure [TS] demonstrates 
the predicted radial profile of the anisotropy parame- 
ter, /3, of the velocity dispersions first and second star 
particles as compared to the overall distribution. Here, 
0=1 — 0.5of/of, where ay is the radial velocity dis- 
persion and of = o-g + a^ is the tangential velocity dis- 
persion. The bin sizes are adjusted so that each consists 
of 100 particles for first and second star case and 1000 
particles for the overall distribution. Here, we focus on 
the results within 50 kpc. Overall DM particles show 
almost isotropic velocity dispersion in the central region, 
and (3 gradually increases, i.e., the radial velocity dis- 
persion becomes more dominant, with increasing radius. 
This is consistent with previous equally high-resolution 
N-body simulations of a Milky Way size halo (e.g., Dic- 
mand et al. 2005). Compared with the overall particles, 
the first star particles have higher /3, i.e., the radial ve- 
locity dispersion is more dominant, especially at smaller 




Fig. 14. — The eccentricity, e = ( f'max — ?*min)/ (^max "t" ^min); as a function of the radius at z — for first (left panel) and second (center 
panel) star praticles, contrasted with the full distribution (right panel). In all cases the line gives the average value as a function of radius. 



radii. This is expected, because the first star particles 
in the inner region form in halos closer to the density 
peak of the main system at high redshifts, and fall into 
the main system at earlier epochs. Second stars have a 
degree of anisotropy between the first particles and the 
overall distribution. 

We also measure the eccentricity of the orbits of the 
first and second stars in Figure [TU To compute these 
values, we run an N-body simulation for these particles 
in the fixed potential of the final system. Here, we fix 
the position of all the particles within the virial radius, 
and follow the evolution of the first and second particles 
from their position and velocity at z = 0. We then com- 
pute the motion of each particle over a full orbit after its 
first passing through pericenter or apocenter, and mea- 
sure the maximum, r max , and minimum, r m i„, radii over 
this full trajectory. The eccentricity is then defined as, 
e = (r max - r min )/(r max + r min ). Figure [H shows that 
both first and second particles have a relatively high ec- 
centricity ((e) « 0.7). This slightly exceeds the mean 
value of the overall distribution, which is roughly consis- 
tent with the observed value of (e) ss 0.6 (Chiba & Beers 
2000). 

4. CONCLUSIONS 

Theoretical models of PopIII star formation have fo- 
cused on their overall cosmological distribution, or their 
contribution to the global star formation history of the 
Galaxy. But z — observations are not of this type. 
As PopIII stars would have been carrying out nuclear 
burning for roughly a Hubble time, they must be faint, 
low-mass stars that are only observable nearby, and in en- 
vironments without significant dust extinction or crowd- 
ing from higher-mass stars. For this reason searches for 
metal-poor stars have been targeted only in the halo of 
the Milky Way, with an emphasis on halo stars roughly 
in the solar neighborhood. This means that comparisons 
between theory and observations have been dependent 
on extremely uncertain extrapolations, drawn from mea- 
surements in a limited environment. 

In this work, we have carried out the first theoretical 
investigation able to quantify the impact of this extrap- 
olation. Our method relies upon two key ingredients. 
The first of these is an extremely high-resolution N-body 
simulation of the Milky Way, which is capable of reliably 
capturing the formation of objects down to the 10 4 K 
atomic cooling limit. The second is a simple and flex- 



ible analytic model of metal-dispersal, which allows us 
to quickly identify regions of metal-free star formation, 
despite the enormous uncertainties in the properties of 
the first stars and supernovae. 

Our results are robust and optimistic. Despite chang- 
ing the efficiency of metal dispersal by two orders of 
magnitude, adopting widely different estimates for the 
distances between neighboring objects, and varying our 
assumptions as to the regions of star-formation within 
collapsed objects, local measurements place substantial 
constraints on PopIII star formation in all of our models. 
In particular we have found that: 

(1) Galactic PopIII star formation occurs over a 
range of redshifts, as metal pollution takes a substan- 
tial amount of time to affect all the Milky Way progen- 
itors. While Galactic PopIII star formation is likely to 
have peaked at z k, 10, it should have continued at ap- 
preciable values down to z k, 5, paralleling the overall 
cosmological evolution. 

(2) The mass fraction of PopIII stars in the Galactic 
halo near the solar orbit is only slightly lower than it is in 
the bulge. Thus, if they have sufficiently long lifetimes, 
a significant number of stars formed in initially primor- 
dial star clusters should be found in the Galactic halo. 
This means that there is no compelling theoretical rea- 
son to motivate observational searches in more difficult 
environments, and present observations should be taken 
as directly constraining the properties of PopIII stars. 
This picture is also consistent with observational studies 
of the lithium-7 abundance of metal-poor stars, which 
suggests that a large fraction of the baryonic matter in 
the early halo may have been proceed by PopIII stars. 

(3) While our models are not able to quantify self- 
enrichment within individual star clusters in detail, we 
expect that stars enriched purely by metal-free stars in 
the same cluster will exist cospatially with the "first 
stars" in our model. Furthermore we are able to make 
statements about the distribution of "second stars" that 
do not form cospatially with PopIII stars, but rather 
form out of gas that is enriched by neighboring clusters 
of PopIII stars. In almost all of our models, the mass 
fraction of such stars in the Galactic halo is compara- 
ble to that of PopIII stars, although there are very few 
second stars in our weakest £^ n — 10~ 4 , Eulerian case. 
Unless the efficiency of PopIII metal ejection was at the 
lower end of the allowed range, many of these stars should 
also be found nearby. 
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(4) We find that the final distribution of first and sec- 
ond stars remains practically unchanged when we move 
from our fiducial model, which associates stars with all 
dark-matter particles in each first and second generation 
object, to a centrally-concentrated model, which asso- 
ciates stars with only the inner 10% of particles in each 
object. Thus we do not expect the spatial distribution 
of stars within each first and second generation object to 
affect their distribution today. 

(5) Incorporating the effects of reionization leads to 
an overall smaller number of satellite galaxies containing 
first and second stars, and helps to address the "missing 
satellite problem." However including such effects does 
not change the overall widespread distribution of these 
stars. 

While stars are the sources of all metals, metals need 
not be contained in every star. Furthermore, we have 
shown here that stars without metals need not to be lo- 
cated near the Galactic center. Rather observations of 
metal-poor halo stars are providing us with direct con- 
straints on the properties of the first stars. Though pri- 
mordial stars were formed in the distant past, some of our 



greatest clues to the process of primordial star formation 
are likely to come from our local Galactic neighborhood. 
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